'''
Created on 22.8.2012

@author: mphkh
'''
from scipy.stats import norm
from math import log
from math import sqrt
from math import exp


def bs_d1(Price, Strike, Rate, Time, Volatility, Yield):
    return (log(Price/Strike)+(Rate - Yield + Volatility**2/2)*Time)/(Volatility*sqrt(Time))

def bs_d2(Price, Strike, Rate, Time, Volatility, Yield):
    return float((log(Price/Strike)+(Rate - Yield + Volatility**2/2)*Time)/(Volatility*sqrt(Time)) - Volatility*sqrt(Time))

def blackscholes(Price, Strike, Rate, Time, Volatility, Yield):
    d1 = bs_d1(Price, Strike, Rate, Time, Volatility, Yield)
    d2 = d1 - Volatility*sqrt(Time)
    
    cprice  = Price*exp(-Yield*Time)*norm.cdf(d1)-Strike*exp(-Rate*Time)*norm.cdf(d2)
    pprice = Strike*exp(-Rate*Time)*norm.cdf(-d2)-Price*exp(-Yield*Time)*norm.cdf(-d1)
    
    return [cprice, pprice]

def blackscholesdelta(Price, Strike, Rate, Time, Volatility, Yield):
    d1 = bs_d1(Price, Strike, Rate, Time, Volatility, Yield)

    cdelta = delta = exp(-Yield*Time)*norm.cdf(d1)
    pdelta = exp(-Yield*Time)*(norm.cdf(d1)-1)

    return [cdelta, pdelta]

def blackscholesgamma(Price, Strike, Rate, Time, Volatility, Yield):
    d1 = bs_d1(Price, Strike, Rate, Time, Volatility, Yield)

    return norm.cdf(d1)*exp(-Yield*Time)/(Price*Volatility*sqrt(Time))

def blackscholesvega(Price, Strike, Rate, Time, Volatility, Yield):
    d1 = bs_d1(Price, Strike, Rate, Time, Volatility, Yield)

    return Price*exp(-Yield*Time)*norm.cdf(d1)*sqrt(Time)



Price = 100.0
Strike = 105.0
Rate = 0.05
Time = 1.0
Volatility = 0.10
Yield = 0.02

a = blackscholes(Price, Strike, Rate, Time, Volatility, Yield)
b = blackscholesdelta(Price, Strike, Rate, Time, Volatility, Yield)
c = blackscholesgamma(Price, Strike, Rate, Time, Volatility, Yield)
d = blackscholesvega(Price, Strike, Rate, Time, Volatility, Yield)

print a[0],a[1]
print b[0],b[1]
print c
print d
